#include "ccd.h"
#include "hfsolver_ho.h"
#include "mbpt.h"
#include "rpa.h"
#include "system_hf.h"
#include <fstream>
#include <iostream>
// #include"adc2.h"
//#include"spsolver.h"
//#include"system_gammow.h"
//#include"sm.h"

using std::cin;
using std::cout;
using std::endl;
using std::ofstream;
using std::to_string;
// using namespace SM;

int main() {
  // System_Tz_HO
  // h("../../../../public/n3lo_vlowk/sp3.dat","../../../../public/n3lo_vlowk/v3.dat");
  // System_Tz_HO h("./sp.dat","./vint.dat");
  //  System_Tz_HO h("./sppn1.dat","./vpn1.dat");
  // System_Tz_HO h("./spox16.dat","./vint2.dat");
  //  System_Tz_HO h("./sp_n.dat","./vint_n.dat");
  //  System_Tz_HO
  //  h("/home/wuqiang/project/vlowk/sppn.dat","/home/wuqiang/project/vlowk/vpn.dat");
  //  System_Tz_HO h
  //  ("/home/public/v2b_n3lo_vlowk/2.1/18/6/06vorbits.dat","/home/public/v2b_n3lo_vlowk/2.1/18/6/06v18.dat");
  //  h.setup();
  // System_Tz_HO
  // h("/home/public/heff_v2b_av18_vlowk/1.9/18/8/08vorbits.dat","/home/public/heff_v2b_av18_vlowk/1.9/18/8/08v18.dat");
  string prefix;
  string orbfile;
  string intfile;
  //  //    cin>>prefix>>orbfile>>intfile;
  /*prefix = string("/home/public/n3lo_vlowk/");
  ////  prefix=string("/home/public/v2b_n3lo_vlowk/2.1/22/7/");
  cin >> orbfile >> intfile;
  ////  orbfile=string("07vorbits.dat");
  ////  orbfile=string("sp_N8_hw22_2.1.dat");
  ////  intfile=string("v_N8_hw22_2.1.dat");
  //  //  intfile=string("07v22.dat");
  orbfile = prefix + orbfile;
  intfile = prefix + intfile;
  string nuclei;
  cin >> nuclei;
  int Nshell;
  cin >> Nshell;
  float lambda;
  cin >> lambda;
  char str_lambda[5];
  sprintf(str_lambda, "%3.1f", lambda);
  string file = nuclei + string("_response_") + string(str_lambda) +
                string("_") + to_string(Nshell) + string(".txt");
  ofstream fout(file.c_str());

  System_Tz_HO h(orbfile.c_str(), intfile.c_str());
  h.setup();
  //  h.printOrbitals();

  // // System_T_HO ht("10jorbits.dat","10j08.dat",8);
  // // ht.setup();

  // //TtoTz(ht,h);

  int A = 4;
  int Z = 2;

  cin >> A >> Z;

  cout << orbfile << endl;
  cout << intfile << endl;
  cout << Nshell << endl;
  cout << lambda << endl;
  cout << A << "\t" << Z << endl;
  cout << "Out put file : " << file << endl;*/

  // SHFSolver_HO hf(&h,A,Z);
  // hf.iter();
  // hf.calEnergy();
  // cout<<hf.Energy<<endl;
  // hf.printOrbitals(10);

  // System_Tz_WS hws("cspox16.dat","cvint2.dat");
  // hws.setup();
  // SHFSolver_WS hfws(&hws,16,8);
  // hfws.iter();
  // hfws.calEnergy();
  // cout<<hfws.Energy<<endl;
  // hfws.printOrbitals(20);

  // System_Tz_K hk(64);
  // hk.setup();
  // SHFSolver_K hfk(&hk,16,8);
  // hfk.iter();
  // hfk.calEnergy();
  // cout<<hfk.Energy<<endl;

  // cout<<"A\tnuclei\tHF\tRPA\tMBPT2\tMBPT3(4p2h)\tMBPT3(2p4h)\tMBPT3(3p3h)\tCCD_Ladder\tCCD\n";
  // cout<<A<<"\t";
  // cout<<nuclei<<"\t";

  /*System_Tz_HF hhf(&h, A, Z); // A,Z
  hhf.setup();*/

  // cout<<"rms:\t"<<hhf.rms()<<endl;
  // fout<<"rms:\t"<<hhf.rms()<<endl;

  /*hhf.printOrbitals();*/

  // cout<<hhf.HFenergy<<"\t";
  // cout<<hhf.radialOverlap(4,5)<<endl;

  ////  for nnlo_sat,change get2B in system_hf.h and get1B in system_hf.cpp
  int Flag_TBME = 2; // ! changed by zhuo
                     // 0: get2BmatOnNeed
                     // 1: setupTwoBodyMat()&get2BmatAt
                     // 2: get2BmatAt only, For NNlosat
  int A = 24;
  int Z = 16;
  System_Tz_HO h;
  System_Tz_HF hhf(&h, A, Z, Flag_TBME);
  string path = string("/home/public/nnlosat/O24/");
  string suffix = string("nnlosat_N12E12_hw22_O24A24.dat");
  string spfile = path + string("sp_energy_") + suffix;
  string umatfile = path + string("umat_") + suffix;
  string kinfile = path + string("kinetic_") + suffix;
  string vfile = path + suffix;
  setup_nnlo_sat(h, hhf, spfile.c_str(), umatfile.c_str(), kinfile.c_str(),
                 vfile.c_str());
  string data_path = "data/O24_nnlosat/";
  // ofstream fout(file2.c_str());
  // cout << nnlo_sat_getHFenergy(hhf) << "\t";
  // MBPT mbpt(hhf);
  // complexd E2p2h,E4p2h,E2p4h,E3p3h;
  // E2p2h=mbpt.diagram0B_2p2h(0);
  // E4p2h=mbpt.diagram0B_4p2h(0);
  // E2p4h=mbpt.diagram0B_2p4h(0);
  // E3p3h=mbpt.diagram0B_3p3h(0);
  // cout<<(E2p2h+E4p2h+E2p4h+E3p3h).real()<<"\t";

  RPA rpa(&hhf);
  // rpa.setTDAflag();
  rpa.cal();

  int type = -1, v = 0;
  // cout<<"F\n";
  // rpa.set_q_for_Beta_F(type);//type
  // cout<<rpa.B_F(type,0)<<endl;//type,v
  // rpa.print_XY_ce(0,0,type,v);//J,Par,Tz,v
  // rpa.print_detail_info_ce(0,0,type);//J,Par,Tz

  /*
    cout << "GT\n";
    rpa.set_q_for_Beta_GT(type);
    cout << rpa.B_GT(type, v) << endl;
    rpa.print_XY_ce(1, 0, type, v);       // J,Par,Tz,v
    rpa.print_detail_info_ce(1, 0, type); // J,Par,Tz

    cout << "sum rule\n";
    cout << rpa.S_F() << endl;
    cout << rpa.S_GT() << endl;
    */

  // cout<<rpa.getDeltaE()<<"\t";
  // rpa.calOccuNum();
  // for(int i=0;i<10;i++)
  //   {
  //     cout<<i<<"\t"<<rpa.OccuNum[i]<<endl;
  //   }
  // fout<<rpa.S(0,0,0)<<endl;
  // cout<<rpa.S(0,0,0)<<endl;
  // cout<<rpa.S(1,1,1)<<endl;
  // cout<<rpa.S(2,0,0)<<endl;
  //  //rpa.print_rou_tr(0,0,12,100);

  rpa.print_detail_info(1, -1, 1); // J,Par,Tz

  string filename_ISM = data_path + "ISM.txt";
  string filename_IVD = data_path + "IVD.txt";
  string filename_ISQ = data_path + "ISQ.txt";

  ofstream file_ISM(filename_ISM.c_str());
  file_ISM << "isoscalar monopole\n";
  rpa.R(0, 0, 0, file_ISM);
  file_ISM.close();
  ofstream file_IVD(filename_IVD.c_str());
  file_IVD << "isovector dipole\n";
  rpa.R(1, 1, 1, file_IVD);
  file_IVD.close();
  ofstream file_ISQ(filename_ISQ.c_str());
  file_ISQ << "isoscalar quadrapole\n";
  rpa.R(2, 0, 0, file_ISQ);
  file_ISQ.close();

  //  rpa.R(1,1,0,cout);
  //  cout<<endl<<endl<<endl;
  // rpa.print_detail_info(1,1,0);
  // rpa.print_detail_info(0,0,0);
  //  cout<<endl<<endl;
  //  rpa.print_detail_info(1,1,1);
  //  cout<<endl<<endl;
  //  rpa.print_detail_info(2,0,0);

  // ERPA erpa(&hhf);
  // erpa.init();
  // erpa.iteration();
  // erpa.R(1,1,1,cout);

  // MBPT mbpt(&h,A,Z);
  // mbpt.setup(A);
  // complexd E2p2h,E4p2h,E2p4h,E3p3h;
  // E2p2h=mbpt.diagram0B_2p2h(0);
  // E4p2h=mbpt.diagram0B_4p2h(0);
  // E2p4h=mbpt.diagram0B_2p4h(0);
  // E3p3h=mbpt.diagram0B_3p3h(0);
  // cout<<E2p2h.real()<<"\t"<<E4p2h.real()<<"\t"<<E2p4h.real()<<"\t"<<E3p3h.real()<<"\t";
  //  cout<<(E2p2h+E4p2h+E2p4h+E3p3h).real()<<"\t";
  //  cout<<mbpt.HFSystem.HFenergy+E2p2h+E4p2h+E2p4h+E3p3h<<endl;

  // mbpt.setPSpace(6);
  // mbpt.setPSpace(7);
  // mbpt.setPSpace(8);
  // mbpt.setPSpace(9);
  // mbpt.setPSpace(10);
  // mbpt.setPSpace(11);
  // // mbpt.setPSpace(12);
  // // mbpt.setPSpace(13);
  // // mbpt.setPSpace(14);
  // // mbpt.setPSpace(15);
  // // mbpt.setPSpace(2);
  // // mbpt.setPSpace(3);
  // // mbpt.setPSpace(4);
  // // mbpt.setPSpace(5);

  // //  ofstream fout("temp",std::ios::app);

  // //  mbpt.printP2BStates();
  // //  mbpt.solve0B();
  // //  cout<<"0B:\t"<<mbpt.CoreSumSPE<<endl;
  // //  mbpt.printHFOrbitals(20);
  // //  mbpt.print0BGraphData(fout,0,0.1,600);//start,step,n
  // mbpt.setup(A+1);
  // mbpt.printHFOrbitals(12);
  // mbpt.printPGroups();
  // //  cout<<"1B:\t"<<mbpt.CoreSumSPE<<endl;
  // // cout<<mbpt.HFSystem.HFenergy<<endl;
  // mbpt.solve1B();
  // //  mbpt.print1BGraphData(fout,0,10,0.1,100);//block,start,step,n
  // mbpt.setup(A+2);
  // mbpt.printHFOrbitals(12);
  // mbpt.printP2BChannels();
  // mbpt.solve2B();
  // //  cout<<"2B:\t"<<mbpt.CoreSumSPE<<endl;
  // //  mbpt.print2BGraphData(fout,11,-30,0.1,500);//block,start,step,n
  // mbpt.printResult();
  // mbpt.mat2Bminus1B();
  // //  mbpt.symmetrize();
  // mbpt.printP1BMat();
  // mbpt.printP2BMat();

  //  mbpt.printPOrbitals();
  //  mbpt.printP1BMat();
  //  mbpt.printP2BMat();
  //  cout<<mbpt.CoreSumSPE<<endl;
  // for(int i=0;i<40;i++)
  //   {
  //     double w=0+i;
  //     //
  //     cout<<w<<"\t"<<mbpt.get1BQBox(w).real()<<"\t"<<mbpt.get1BZBox(w).real()<<endl;
  //     MBPT::Mat mat;
  //     mbpt.get2BQbox(0,mat,w);
  //     cout<<w<<"\t"<<mat(0,0).real()<<"\t";
  //     mbpt.get2BZbox(0,mat,w);
  //     cout<<w<<"\t"<<mat(0,0).real()<<"\t"<<endl;
  //   }

  // mbpt.debug();
  // mbpt.calOccuNum();
  // double num=0;
  // for(int i=0;i<mbpt.HFOrbs.size();i++)
  //   {
  //     cout<<mbpt.OccuNum[i]<<endl;
  //   }

  // CCD ccd(&hhf);
  // // cout<<ccd.TwoOrder_2p2h()<<endl;
  // if(ccd.Ladders())
  //   cout<<ccd.CCD_deltaE<<"\t";
  // if(ccd.Full())
  //   cout<<ccd.CCD_deltaE<<"\t";
  // cout<<endl;

  // ADC2 adc(&hhf);
  // adc.setup();
  // adc.solveDysonMatrix_Lanczos(3,3);
  // adc.calEnergy();
  // cout<<adc.Energy<<endl;
  // cout<<adc.num<<endl;
  // adc.printResults();

  //  SPSolver spsolver(&h,4,2);
  //  spsolver.init();
  // //  spsolver.printVrr(2);
  // //  cout<<spsolver.getVkkfromHF(0,0,0)<<endl;
  //  spsolver.solve();
  // cout<<spsolver.Energy<<endl;
  //  spsolver.printOrbitals();
  //  spsolver.setupOverlap(1);
  //  spsolver.printOrbitals(-1,5);
  // spsolver.printOrbitals(-1,-1);//num,channel, -1 print all
  // for(int i=0;i<100;i++)
  //   {
  //     double k=0.1*i;
  //     cout<<k<<"\t"<<spsolver.hf_k(1,k).real()<<endl;
  //   }
  //  spsolver.printVkk(0);

  // System_Tz_Gammow gam(&h,4,2);
  // gam.setup();
  // gam.printOrbitals();
  //  gam.print1B("1B.dat");
  //  gam.print2B("2B.dat");

  // System_Tz_SM effV("usdbpn.sp","usdbpn.int");
  // effV.setup();

  // System_Tz_SM_M effVm(&effV);
  // effVm.setup();
  //  effVm.printOrbitals();

  // Configs configs(&effVm,2,0);//,vN,vZ
  // configs.setup(0,0);//MM,Par
  // cout<<configs.getME(0,0)<<endl;
  // cout<<configs.getJ2(1,1)<<endl;

  // SMSolver smsolver(&effVm,2,0);//,vN,vZ
  // smsolver.diag(0,0);//MM,Par
  // smsolver.printStates(10);
  // smsolver.diag_lanczos(10,0,0);//num,MM,Par,step=1
  // smsolver.printStates(10);
  // smsolver.calOccuNums(0);
  // smsolver.printOccuNums();
  return 0;
}
